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ABSTRACT 

The spatial-frequency coverage of a radio interferometer is increased by combining samples acquired at different 
times and observing frequencies. However, astrophysical sources often contain complicated spatial structure that 
varies within the time-range of an observation, or the bandwidth of the receiver being used, or both. Image 
reconstruction algorithms can been designed to model time and frequency variability in addition to the average 
intensity distribution, and provide an improvement over traditional methods that ignore all variability. This 
paper describes an algorithm designed for such structures, and evaluates it in the context of reconstructing 
three-dimensional time-varying structures in the solar corona from radio interferometric measurements between 
5 GHz and 15 GHz using existing telescopes such as the EVLA and at angular resolutions better than that 
allowed by traditional multi-frequency analysis algorithms. 
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1. INTRODUCTION 

Image reconstruction in radio interferometry^ benefits from maximizing the number of distinct spatial frequencies 
at which the visibility function of the target is sampled. For a given array, the spatial-frequency coverage can 
be increased by taking measurements across several hours (Earth-rotation-synthesis) and at multiple observing 
frequencies (multi- frequency-synthesis). However, astrophysical sources sometimes show time and frequency 
variability within the parameters of such an observation. Reconstruction algorithms have traditionally ignored 
this variability, leading either to imaging-artifacts or to data-analysis strategies that use subsets of the data 
independently, both which limit imaging sensitivity and accuracy. 

In recent years, broad-band receivers have been installed on radio interferometers specifically to increase 
imaging sensitivity, and this has necessitated the development of algorithmJ^Sl that account for sky spectra in 
addition to average intensity. A similar idea applies to time- variability as well, and recent worlPhas demonstrated 
combined time and frequency modeling for point-sources (without multi-scale support), mainly for the purpose 
of eliminating artifacts from images of the average intensity distribution. 

This paper presents a reworking of the system of equations being solved when modeling time and frequency 
varying structure, and includes multi-scale support. Using a software implementation of a multi-scale, multi- 
frequency, time-variable image-reconstruction algorithm that can be operated in modes that emulate various 
intermediate algorithms, a series of imaging tests were carried out on data simulated for three-dimensional time- 
varying structure in the solar corona. The results are discussed in terms of the accuracy of reconstructing the 
time and frequency-dependent structure at an angular resolution better than what traditional methods allow. 

2. MODEL-FITTING TO RECONSTRUCT AN IMAGE 

According to the van-Cittert-Zernicke (VCZ) theorem^ , an ideal radio interferometer samples the visibility func- 
tion of the sky brightness distribution. For a single frequency channel v and integration timestep t, each pair of 
detectors corresponds to one visibility sample, located at the spatial-frequency defined by the distance between 
the pair, the phase-tracking center, and the observing frequency. With n ant antennas, n ant (n an t l)/2 such mea- 
surements of the visibility function are made at spatial frequencies that define the uu-coverage of the instrument. 
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Image-reconstruction is then a non-linear process that solves for the parameters of a sky model, subject to con- 
straints from the data measured at this set of spatial frequencies. The instantaneous single-channel ww-coverage 
is usually augmented by Earth-rotation synthesis and multi-frequency synthesis. Under the assumption that 
the sky structure is invariant with time and frequency, these samples can be considered additional indepen- 
dent measurements of the same visibility function. The VCZ theorem holds, and standard image-reconstruction 
algorithms such as CLEAN 6 3 and MEMEE apply, and benefit from the extra data-constraints. 

When the sky structure varies with frequency or time, within the parameters of a single observation, the 
traditional solution has been to partition the data in time and frequency, such that the VCZ theorem holds on 
each piece. This approach sacrifices sensitivity and ui>-coverage during the non-linear reconstruction process. 
While this suffices for simple point- like structures and emission with high signal-to-noise ratios, non-uniqueness 
in the reconstruction of complicated emission can introduce errors that are inconsistent across time and/or 
frequency. In such cases, it is better to model and fit for the frequency and time-dependence at the same time 
as the average intensityP Although this approach has primarily been applied to separate the varying and non- 
varying components of the dataset with the goal of producing a single image representing non- variant structure, 
it is sometimes possible to reconstruct the time and/or frequency variation with enough accuracy to be used for 
astrophysical interpretation. Such an analysis is strongly dependent on the models used during reconstruction, 
how well the measurements constrain their parameters, and how well the reconstruction algorithm is able to 
converge to the best-fit solution. 

For several decades, the general approach to image reconstruction in radio interferometry has been to describe 
the sky structure as a collection of delta-functions, parameterized by their location and amplitude.^ For sky 
emission comprised of several discrete sources, delta-functions are a compact representation of the signal being 
modeled. More recently, multi-scale algorithms^HSl have been used to model extended emission, by decomposing 
the sky brightness into a linear combination of basis functions representing flux-components of different shapes 
and sizes. This is considered a sparse representation of extended spatial structure, and the parameters to be 
solved-for are the location, amplitude, and size of the components. Multi-frequency algorithms 2 4 model and 
reconstruct the sky spectrum by solving for the coefficients of a Taylor-polynomial representation of the spectrum, 
again, a sparse representation of smooth spectral structure. This idea can be extended to time-variability^ as 
well, in which a set of basis functions is chosen so as to compactly represent the time-variation. 

In most of these cases, the image-reconstruction algorithm first transforms the data into a basis in which the 
model has a sparse representation, and then applies a greedy algorithm to search for flux components and their 
best-fit parameters. 

MS-MF- TV-CLEAN, an algorithm that uses a sparse representation of multi-scale, frequency-dependent, 
time-variable structure and solves for its parameters via x 2 minimization is described below. Differences with 
the multi-frequency and time- variable algorithms in Refs. 3, 5} |10| and the consequences of these differences will 
be pointed out when relevant. Sec. [3] uses some imaging results to evaluate the accuracy at which such methods 
are able to reconstruct time and frequency-variable structure, using uv-coverages offered by moderately sized 
interferometers such as the EVLA with 27 antennas and a scaled-down FASR with 15 antennas. 

2.1 Multi-term image models 

Parameterized models for multi-scale, frequency-dependent and time-varying spatial structure, and their combi- 
nation as used in the MS-MF- TV-CLEAN algorithm are described below. 

Multi-Scale structure : For a finite set of N s spatial scales, a multi-scale image model is written as 

N a -1 

T m = if p *if v (i) 

s=0 

where I" hp are multi-scale basis functions, and I s s ky is an image of delta-functions that mark the locations and 
amplitudes of the basis function of type s. As implemented in MS-CLEAN 10 , the first scale function J*_g is 
chosen to be a 5-function in order to always allow for the modeling of unresolved sources, and successive basis 
functions are inverted and truncated parabolas of increasing width (as s increases) . 



Multi-Frequency structure A sky brightness distribution that varies smoothly with observing frequency can 
be modeled with Taylor-polynomials per pixel, and represented as a linear combination of coefficient images and 
basis functionsP^ The intensity at each frequency channel can be written as 

where N p is the total number of terms in the series and I p ky ; < p < N p are the coefficient images (one set of 
coefficients per pixel), v is the observing frequency, v$ is a chosen reference frequency. 



Time- Variable structure A smooth variation of the source intensity as a function of time can be described 
as another Taylor-polynomial. 
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where N q is the total number of terms in the series and I q v ; < q < N q are the coefficient images (one set of 
coefficients per pixel), t is the time, and t m id,t max ,t m i n are calculated from the time-range of the observation. 
These basis functions correspond to a Taylor polynomial evaluated between -1 and +1. Ref. |5| discusses and 
demonstrates that if the time variation is not smooth, a Fourier-basis is much more appropriate, but requires 
many more terms in the series. 



Multi-scale multi-frequency time-variable structure The basis functions described in Eqns. (??) and Q 
can be combined as follows. 
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where the expression inside the square brackets represents images of delta functions that mark the locations of 
flux components of shape I* hp and whose amplitudes follow dynamic-spectra defined by linear combinations of 
polynomials along the time and frequency axes. 

With N s multi-scale basis functions, N p Taylor-terms along frequency, and N q Taylor-terms along time, the 
total number of basis functions becomes N s x (N p + N q — 1). Alternately, one could use Zernicke polynomials to 
describe patterns on the 2D time-frequency plane. 

Note that for each scale-size, this choice of basis-set requires fewer terms than the N p x N q basis functions 
used in the multi-beam clean algorithm described in Ref. [5] which also uses a Gram-Schmidt orthogonalization to 
regularize the problem. Also, the algorithm described here includes a multi-scale model, which from the analysis 
done to develop the MS-MF-CLEAINP becomes very relevant when attempting to use all the reconstruction 
products for astrophysical analysis. 

2.2 Multi-term measurement equation 

Eqn. @ can be combined with the measurement process of an imaging interferometer to give the following 
Measurement equations (written in matrix notation). 
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where I r is an m x 1 list of pixel amplitudes that represents {1^ ky , I s s k p v , T s s k q v Vs,p, q} as written in Eqn. Q. 
Each [A r ] is an n x m measurement matrix for the quantity I r which produces a list of n visibility measurements 



(n = nbaseiines x n t i m esteps x Tichannels)' Ignoring direction-dependent instrumental effects, we can write [A r ] = 
[W r ] [S] [F] where [F] is an m x m Fourier-transform operator, [S] is an n x m matrix of ones and zeros representing 
the uv-coverage for all baselines over frequency and time, and [W r ] is an n x n diagonal matrix of weights that 
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The measurement equations in block matrix form (for N r — 3) are given by 



[A ] [A,] [A 2 



lo 
h 
h 



V 



(6) 



All N r measurement matrices of shape n x m are placed side-by-side to form a larger measurement matrix. The 
list of parameters becomes a vertical stack of N r vectors each of shape m x 1. The new measurement matrix of 
shape n x mN r operates on an mN r x 1 list of parameters to form annxl list of measurements. 

A least-squares solution of this system of equations requires the normal equations to be constructed and 
solved as described below. 

2.3 Multi-term normal equations 

The normal equations (for N r = 3) are given by 
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where [W] is an n x n diagonal matrix of data weights, usually based on system noise levels. The Hessian matrix 
on the LHS consists of N r x N r blocks, each of size m x m. The list of parameters is a stack of N r column 
vectors I r , and the RHS vector is a set of N r weighted inversions of the data vector. 

When [A r ] — [W r ][S][F], each Hessian block [A' Ti W A r . ] is a Toeplitz convolution operator, with each row 
containing a shifted version of a point-spread-functiorllMiff . Each image vector in the RHS is therefore a sum 
of convolutions of different image pixel vectors and convolution kernels, and the process of reconstructing I r Vr 
is equivalent to a joint deconvolution. 

This system of equations is evaluated by first computing the RHS vectors and one point-spread-function per 
Hessian block. 

The point spread function for Hessian block ri,rj is computed as [F^S^W rt WW r Al- This translates to the 
inverse Fourier transform of the product of the uw-coverage, data weights [W] and the spatial-frequency signature 
of the basis functions labeled by , rj . 

The RHS vector contains a stack of N r residual images, each computed as [•FT J SrWT r W]y'. This translates to 
weighting the visibilities using data weights as well as weights that evaluate the r th basis function (a matched- 
filtering step), resampling the visibilities onto a regular grid and then taking an inverse Fourier transform. It is 
important to note that basis-function weights for the time and frequency axes must be applied in the spatial- 
frequency domain, prior to resampling the list of visibilities onto a regular grid. This is to ensure that if multiple 
visibility values from different baselines, frequencies and times get mapped onto a single uv grid cell, all the 
contributing visibilities have the correct weights applied to them. Weights for the multi-scale basis functions, 
however, can be applied as a multiplication on the gridded visibilities. 



Differences with the Sault-Wieringa approach : The Sault-Wieringa multi-beam algorithms of Refs. [3] 
and [5] constructs RHS vectors by convolving the standard residual image with point-spread- functions constructed 
for each series term (I^ hs = lQ hs */p s ^). fn the visibility-domain, this corresponds to the multiplication of the 
gridded visibilities with gridded weights. If there are any uv grid cells with contributions from measurements 
from different frequencies, they do not get their correct individual weights. The Sault-Wieringa approach works 
well with interferometers with minimal overlap in spatial-frequency coverage (as demonstrated by Ref. |3jwith the 
ATCA and Ref. 5 with e-MERLIN, but incurs errors for arrays with dense spatial-frequency coverage such as the 
EVLA. For multi-scale multi-frequency imaging, changing the computations to follow those discussed here and 
in Ref. 4 eliminated this instability and allowed reconstructions of the source spectrum for extended emission at 
an accuracy sufficient for astrophysical analysis. A similar argument holds for the SW-based algorithm described 
in Ref. [5] when solving for time- variability. 

2.4 Solving the Normal equations 

The normal equations shown in Eqn. ^ are solved iteratively. 

First, the pseudo-inverse of the Hessian is computed by inverting a block-diagonal approximation of the 
matrix for the time-dependent and frequency-dependent terms and a diagonal approximation for the multi-scale 
terms (required for numerical stability^) . 

This pseudo-inverse is applied to all pixels of the RHS images. A set of N r solutions is then picked by 
finding the location and scale-size corresponding to the flux-components that produces the largest reduction in 
X 2 in the current step. The residual images are updated by evaluating the LHS for each set of components and 
subtracting from the RHS. Iterations of the above steps are performed until the peak residual of the zero-th 
order point-source image reaches the level of the first sidelobe of the zero-th order point-spread function. Model 
visibilities are calculated by evaluating Eqn. residual visibilities and RHS residual images are constructed, 
and the process repeats until residuals are noise-like. These steps are a direct extension of the normalization, 
peak-finding and update steps of the standard CLEAPp! algorithm, and are very similar to the steps of the 
MS-CLEAlNp] , SW-CLEA1NP and SW-Multi-BeanP algorithms. 



Software Implementation : The MS-MF-TV-CLEAN algorithm described here was implemented as an 
extension of the existing MS-MFS algorithm using the CASA software libraries. The following acronyms will be 
used to distinguish between various subsets of this implementation used for the tests in Sec. [3] 

1. MS-CLEAN : The model contains only multi-scale terms. The resulting algorithm differs from that in 



Ref. 10 in that it does not require a scale-dependent bias-factor. 



2. MS-MF-CLEAN : The model contains frequency-dependent terms with multi-scale components^ . 

3. MS-TV-CLEAN : The model contains time- variable terms with multi-scale components. The time- variability 
part of the algorithm differs from that in Ref.|5]in its choice of basis functions and how it computes residual 
images. 

4. MS-MF-TV-CLEAN : The model contains time and frequency-variable terms with multi-scale components. 
The MF-TV part of this method differs from that in Ref. [5]in the way the time and frequency basis functions 
are combined and uses fewer total terms in the series expansion of the model. 



3. IMAGING 3D STRUCTURE IN THE SOLAR CORONA 

Gyro-synchrotron radio emission from the solar corona within the frequency range of 1 GHz to 20 GHz traces 
magnetic-field strengths at difference heights above the surface of the surP^ . Details about this 3D magnetic- 
field structure and how it varies with time are essential to the understanding of the physics driving solar flares 
and other coronal activity. Radio interferometric observations at cm wavelengths provide crucial information for 
theoretical models, and several recent efforts^H^ni have made progress on forward-modeling using reconstructed 
multi-frequency image cubes. 
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Figure 1. The snapshot narrow-band spatial- frequency coverage at 10 GHz (LEFT) is compared with the broad-band 
coverage between 5 GHz and 15 GHz, sampled at intervals of 50 MHz (RIGHT). This coverage corresponds to a 15 element 
array, following a randomized log-spiral pattern with upto 400m baselines. This narrow-band coverage is considered sparse 
for structures seen for solar flare loops, and any algorithm that uses the broad-band coverage will have an advantage. 



The image-reconstruction itself, however, has been limited to traditional methods of making images at multi- 
ple separate frequencies, using narrow-band data, and using data accumulated over time-durations shorter than 
the typical variation timescale. With such narrow-band, snapshot interferometric data, the non-linear recon- 
struction process has to work with a limited number of samples that is often insufficient for an un-ambiguous 
reconstruction of the complicated extended emission present in active coronal regions. Further, when data from 
different frequencies are treated separately, the intrinsic frequency-dependent angular resolution of the imaging 
instrument comes into play, and all further analysis of the resulting image cubes must be done at the angular 
resolution offered by the lowest frequency. This can be improved using a wideband imaging algorithm^ that 
reconstructs frequency-dependent structure, but it is still restricted to only snapshot ww-coverages. 

An algorithm such as that described in this paper has the potential of reconstructing spatially-extended 
structure that varies smoothly in shape with both time and frequency within the parameters of a synthesis 
observation. Such an algorithm would be able to take advantage of the increased imaging-fidelity that comes 
from the better spatial-frequency coverage generated by Earth-rotation-synthcsis and multi-frequency-synthesis, 
and not be overly restricted by the intrinsic angular resolution at different frequencies. 

The following sections describe two examples in which radio interferometric observations of 3D structures in 
the solar corona are simulated, and then reconstructed using both traditional and new methods. The goal of 
these tests is to assess if the image-models described in Sec. [5] suffice to reconstruct time and frequency-dependent 
extended-structure at an accuracy sufficient for further analysis, at an angular resolution better than that given 
by the lowest measured frequency. 

3.1 Frequency-dependent structure - a solar flare loop 

Interferometric data for a solar-loop model were simulated for a snapshot observation with a 15-element array 
with broad-band receivers covering the range of 5 GHz to 15 GHz. The antennas follow a randomized log-spiral 
pattern with up to 400m baselines and represents the inner region of the proposed FASR telescope^ . Fig. [l] 
compares the spatial-frequency coverage of the interferometer at 10 GHz with the wideband spatial-frequency 
coverage, both for a single snapshot observation. The intrinsic angular resolution of this instrument is controlled 
by the largest sampled spatial-frequency and varies from 40 arcsec at 5 GHz to 10 arcsec at 15 GHz. 

Figure [2] compares imaging results using traditional and new methods, and compares them with the true 
images used in the simulation. 

The left column shows the true model of a solar-flare loop at 6 different frequencies. This model was 
constructed such that frequencies around 6 GHz trace structure near the top of the solar-flare loop, and higher 
frequencies trace structure deeper in the solar corona, down to the 'feet' of the solar-flare loop. 
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Table 1. This table lists off-source RMS values in the reconstructed average intensity images for the five scenarios 
discussed in Sec. 3.2 and illustrated in Fig. [3] Algorithms that model time and/or frequency variant structure are able to 
take advantage of the extra spatial-frequency coverage offered by the combining the data, and reduce the level at which 
imaging-artifacts occured. 



The middle column shows the reconstruction using the MS-CLEAN algorithm separately on each frequency- 
channel. This is the traditional form of analysis of such wide-band interferometric data, and the image fidelity 
is clearly limited by the varying angular resolution of the instrument as a function of frequency. 

The column on the right shows reconstructions built from a multi-scale multi-frequency model whose coeffi- 
cients were fitted using the MS-MF-CLEAN algorithm on the data combined over all frequencies. With a 3-term 
Taylor-polynomial, this reconstruction was able to recover the 3D structure at an angular resolution of 10 arcsec. 
For the top part of the loop, the algorithm chose flux components whose amplitudes decreased with increasing 
frequency, and for the feet, it chose components whose amplitudes were low at lower frequencies and increased 
with frequency. 

3.2 Time and frequency variable structure - an expanding active coronal region 

A model of an active coronal region was generated such that it shows compact structure in deeper layers (probed 
by higher frequencies), with the upper layers containing an extended plume. As a function of time, these upper 
layers expand spatially, but the lower layers remain relatively unchanged. Data were simulated for the EVLA 
between 5 GHz and 15 GHz, to test whether such an experiment would be feasible with an existing telescope. 
A set of 8 snapshots were simulated, spread across 4 hours. 

Figure [3] shows iw-coverages of various subsets of the data, as well as the average- intensity images produced 
by applying different algorithms to it. Table [T] lists the off-source RMS levels with the different methods, along 
with the number of data-points and the algorithm used in the reconstruction. The true model used in the 
simulations is shown in the first two columns of Fig. |4j 

The results can be summarized as follows. The traditional method of imaging each channel and snapshot 
separately is clearly limited in its dynamic-range and fidelity, even if MS-CLEAN is applied. Also, if MS-CLEAN 
is applied to the entire dataset without accounting for time and frequency variations, artifacts still dominate the 
residuals and affect the on-source reconstructionas well. Intermediate methods of solving for time- variant spatial 
structure one channel at a time (MS-TV-CLEAN), or solving for frequency-dependent structure one timestep 
at a time (MS-MF-CLEAN), yield better results, but the limited uu-coverage still results in imaging artifacts. 
Finally, a combined reconstruction using MS- MF- TV-CLEAN produces the lowest residuals, although even these 
numbers are more than an order-of-magnitude larger than theoretical point-source sensitivity estimate for this 
simulated dataset. The dynamic range of this image is a few hundred to one thousand, which is likely to be 
sufficient for strong radio emission from the sun. 

Figure [4] compares the true source structure with that reconstructed via MS-MF-TV-CLEAN with 2-terms 
along the frequency axis and 2-terms along the time axis, and 5 different spatial scales. The images in the 3rd 
and 4th columns show that most of the time- variable 3D structure has been recovered. This is a promising result, 
but the reconstructions at the highest frequencies show some spurious extended emission. This points to either 
the inadequacy of a 2-term model for the source 'spectrum', or the inherent ambiguity^ in the spectral signature 
at spatial scales that are not sampled at all observed frequencies. 



4. DISCUSSION 



Image-reconstruction algorithms for radio intcrfcrometry that model time and frequency structure can be used 
not-only to improve the fidelity of the average image, but to also reconstruct the frequency-dependent and 
time-variable structure at an accuracy sufficient for astrophysical interpretation. Smooth frequency-dependent 
structure can be reconstructed across frequency ranges wider than previously demonstrated, achieving higher 
angular-resolution than that given by the lowest frequency in the measured band. Reconstructing time- variability 
is also possible, but since angular-resolution is not of concern here, this feature may be useful only in the context 
of improving signal-to-noise and image-fidelity by using all measurements together. 



The example in Sec. 3.2 is a proof-of-concept demonstration that such imaging of frequency-dependent and 
time-variable structures in coronal active regions may be possible with current telescopes such as the EVLA. 
Future telescopes such as FASR will have more than enough spatial-frequency coverage for an unambiguous 
reconstruction at each timestep and channel, even after tapering data from all frequencies to a common low 
angular resolution. 

From a practical standpoint, the MS-MF- TV-CLEAN implementation described in this paper improves upon 
the multi-beam algorithm in Ref . [5] in that it does not suffer from instabilities incurred when there is significant 
overlap between sampled spatial-frequencies. It suggests the choice of N p + N q — 1 basis functions instead of the 
N p x N q , and also folds-in multi-scale support in a way that allows a wide choice of spatial basis functions whose 
amplitudes are given by time and frequency polynomials. 

There is always a trade-off between a basis set that optimally represents the structure being reconstructed 
and requires a minimal number of components to model it, and a basis set whose members the instrument can 
optimally distinguish between (an orthogonal basis for the instrument) but may require many more components 
to model the target structure. For the situations explored in this paper, very simple polynomial basis functions 
produce usable reconstructions, although some care must be taken to ensure that the instrument can distinguish 
between them. This information is encoded in the diagonals of each block of the Hessian matrix, and scale-sizes 
and polynomial orders must be chosen such that its condition number is not too large. Fine structures in the 
dynamic spectra can be reconstructed only with higher-order terms which are increasingly harder for the instru- 
ment to distinguish from one another, and as demonstrated in Ref. [5] may require an explicit orthogonalization 
step. However, this orthogonalization must be restricted to only the time and freqency basis functions because 
the multi-scale basis functions are usually chosen based on structures that optimally describe the structure being 
imaged, and a basis function orthogonalization may result in a much larger number of components being required 
to fit the data. 

One problem with this approach is the current inability to derive accurate error-estimates on the fitted 
parameters, because it depends not-only on the basis functions and the instrument (as encoded in the Hessian and 
covariance matrices), but also on the appropriateness of the chosen model for the structure being reconstructed. 
For instance with MS-MF-CLEAI# , errors in the spectral-index data product are as high as a few hundred 
percent when delta-functions are used to model extended emission. There is also an inherent ambiguity in the 
spectral structure at the largest scales, because of the frequency-scaling of the instruments uw-coverage, and in 
this case, the errors arc dominated by the inability of the measurements to constrain the model and may require 
a-priori informatiorHEE] . To be practically usable, better measures of the reconstruction quality are required. 

The ideas presented in this paper are one step closer to the direct modeling goals presented in Ref. [19] for the 
reconstruction of 3D magnetic-field structures in the solar corona . The next step for such algorithms ought to be 



to step away from pattern-modeling, and fold physical models directly into the reconstruction process. Refs. 18 
and [20] , among others, demonstrate that physical models fitted to single-channel images from an idealized 
interferometer are capable of recovering magnetic-field structures and other information required to undestand 
the physical processes. An algorithm that combines these ideas and reconstructs the physical model directly 
from visibilities would eliminate the intermediate non-linear step of imaging which has its own uncertainties. 
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Figure 2. The frequency-dependent structure of a solar flare loop at 5 GHz, 7 GHz, 9 GHz, 11 GHz and 13 GHz (ROWS 
: Top to Bottom) are shown for the simulated model (LEFT column), reconstructions performed one channel at a time 
(MIDDLE column), and the wide-band reconstruction using the MS-MF-CLEAN algorithm (RIGHT column). Moving 
from lower to higher frequencies (top to bottom) this structure traces the top of the loop and moves deeper in the solar 
corona until the 'feet' of the loop. The angular resolutions of the images in the MIDDLE column range from 52"x28" at 5 
GHz to 14"x7" at 15 GHz, resulting in a plausible reconstruction only at the higher frequencies. For the RIGHT column, 
the angular resolution is 13" x6" at all frequencies shown, and the structure is recovered at the same fidelity at all of these 
frequencies. These reconstructions used the snapshot uv-coverage shown in Fig.[T]and demonstrate the advantage of using 
the broad-band uv-coverage along with an algorithm that accounts for the changing structure as a function of frequency. 
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Figure 3. This figure compares the results of five reconstructions using different subsets of the data and different 
algorithms. All images are displayed with the same gray-scale/range. Two plots are shown for each example; the 
uv-coverage and the reconstructed average intensity (shown below the uv-coverage plots). The first two examples (top 
two rows) show results of the MS-CLEAN algorithm, ignoring time and frequency variable structure. The first example 
(top left) shows the results with narrow-band snapshot uv-coverage, and represents the imaging fidelity achieved via 
the traditional method of handling time-frequency-variable structure by simply partitioning the data and imaging each 
frequency and timestep independently. The second example (top right) shows the results with all the data, taking 
advantage of Earth-rotation and multi-frequency synthesis but ignoring all variations in time and frequency. The next 
three examples (bottom row) show the results of various Multi-Term algorithms. The first example (bottom left) shows 
an MS-TV-CLEAN reconstruction on one-channel of data while fitting for the time-variability. The second example 
(bottom middle) shows an MS-MF-CLEAN reconstruction on a multi-channel snapshot, where the frequency-structure 
is accounted-for. These two examples show some improvement over pure MS-CLEAN, but still show residuals at a level 
higher than what the dataset allows. The last example (bottom right) shows the MS-MF-TV-CLEAN reconstruction 
where data from all channels and timesteps were used and the algorithm fitted simultaneously for time and frequency 
variations in structure. The noise-level achieved is lower than other examples, but still two orders of magnitude higher 
than the theoretical estimate. Table [T] lists off-source RMS noise levels for all 5 examples. 
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Figure 4. This figure shows the multi-frequency and time-variable structure of the simulated radio emission from a 
sunspot, and compares the true structure (COLUMNS 1,2) with that reconstructed via the MS-MF-TV-CLEAN algorithm 
(COLUMNS 3,4). These results are from the same run represented in last pair of plots in Fig. [4] and used data from 20 
channels and 8 timesteps. In each pair of columns, frequency ranges from 6 GHz (TOP) to 14 GHz (BOTTOM), and 
the time difference between each pair of columns is 4 hours. This structure represents the 3D structure of a sunspot 
whose upper layers are expanding spatially as time progresses. The reconstructions (COLUMNS 3,4) show that most of 
the structure has been reconstructed. The sky model solved-for comprised of a linear combination of multi-scale flux- 
components, each of whose amplitudes had an average intensity term, a gradient along frequency, and a gradient along 
time. This simple model resulted in good reconstructions in most of the frequency-range, but showed errors in modeling 
extended emission at frequencies lower than 6 GHz and higher than 12 GHz. 
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